
########################################################################################
########################################################################################
########                  
######## Peter Bucchianeri
########
######## Replication Code for: 'Is Running Enough? Reconsidering the Conventional Wisdom about Women Candidates'
########
######## Last Updated: May 2017
########
########################################################################################
########################################################################################

rm(list=ls())
options(stringsAsFactors = FALSE)

####
# Set Working Directory
####

wd <- "Enter your working directory here" 
setwd(wd)

####
# Install the following R packages if not already installed
####

# install.packages("ggplot2") 
# install.packages("doBy")
# install.packages("gridExtra")
# install.packages("mvtnorm")
# install.packages("stargazer")
# install.packages("xtable")
# install.packages("plotrix")
# install.packages("TeachingDemos")  
# install.packages("rdd")
# install.packages("interplot")
# install.packages("coin")

library(ggplot2) 
library(doBy)
library(gridExtra)
library(mvtnorm)
library(stargazer)
library(xtable)
library(plotrix)
library(TeachingDemos)  
library(rdd)
library(interplot)
library(coin)


#### Create Folder for Figures if Needed
figure_dir <- paste0(wd, "/Figures/")
if(dir.exists(figure_dir) == FALSE){
  dir.create(figure_dir)
}


#####################################################################################################################
######## US House Data - 1972 to 2010
#####################################################################################################################

dat <- read.csv("Bucchianeri_PB_Replication_Data.csv")
dat$party <- factor(dat$party, levels = c("R", "D"))


#############################
#### Function to plot Higher-Order Polynomial RDDs with Binned Averages
#############################

BinPlot <- function(Data, RegObj, binwidth, runningVar, bandwidth, OutcomeCol, xlab, ylab, title, color, xlim, ylim) {
  BinCenters <- seq(-bandwidth + .5*binwidth, bandwidth - .5*binwidth, by=binwidth)
  BinPoints <- NULL
  BinSDs <- NULL
  Seq <- seq(-bandwidth, bandwidth - binwidth, by=binwidth)
  for(i in Seq) {
    Avg <- mean(Data[runningVar >= i & runningVar < i + binwidth, OutcomeCol], na.rm=T)
    if(nrow(Data[runningVar >= i & runningVar < i + binwidth, ]) > 1){
      SD <- sd(Data[runningVar >= i & runningVar < i + binwidth, OutcomeCol], na.rm=T)
    } else {
      SD <- 0
    }
    BinPoints <- append(BinPoints, Avg)
    BinSDs <- append(BinSDs, SD)
  }
  plot(runningVar, Data[,OutcomeCol], cex=.4, las=1, col="gray50", xlab=xlab, ylab=ylab, xlim=xlim, ylim=ylim, pch=16, main=title)
  abline(v=0, lty=2)
  points(BinCenters, BinPoints, col=color, cex=1, pch=20)
}


###############################################################
##### Local Linear Regression Function -- Returns Beta1 for each bandwidth in Seq
###############################################################

LocalLin <- function(Seq, Sub, Formula, PredProb = F, CI = 1.96, use_alt = FALSE, alt_Form = NA, alt_bw = NA){
  dat_sub <- dat[Sub,]
  Est <- rep(NA, length(Seq))
  High <- rep(NA, length(Seq))
  Low <- rep(NA, length(Seq))
  WinDiff <- matrix(NA, ncol=6, nrow=length(Seq))
  num <- 1
  for(i in Seq){
    this_Form <- Formula 
    ## Adjustment to use diff formula if needed (only relevant for Pre/Post analysis with insufficient incumbency variation)
    if( use_alt == TRUE ){
      if( i <= alt_bw ){
        this_Form <- alt_Form
      }
    }
    iSub <- dat_sub[abs(dat_sub$PrimMargin) <= i/2,]   #i/2: 1% BW =  -1/2% to 1/2% // Can do just i for 1% BW = -1% to 1%
    Local <- lm(this_Form, data = iSub)
    Est[num] <- coef(Local)[2]
    RegSE <- summary(Local)$coefficients[2, 2] #Standard SEs
    Robust <- coeftest(Local, vcov = hccm(Local, type="hc1"))[2,2] #Robust SE's
    Low[num] <- Est[num] - CI*max(RegSE, Robust) #Using maximum of regular and Robust SE's
    High[num] <-Est[num] + CI*max(RegSE, Robust)
    ### Predicted Probabilities
    if(PredProb == T){
      Sim <- rmvnorm(10000, mean=coef(Local), sigma=vcov(Local))
      WinPre <- Sim[, "primary_winner"] #Female candidate Wins Primary Pre 1990
      WinPost <- Sim[, "primary_winner"] + Sim[, "Post90"] + Sim[, "primary_winner:Post90"]  #Female candidate Wins Primary Post 1990
      LosePost <- Sim[, "Post90"] #Male Candidate Wins After 1990
      WinDiff_i <- WinPost - WinPre
      WinDiff[num, 1:3] <- cbind(mean(WinDiff_i), quantile(WinDiff_i, .025), quantile(WinDiff_i, .975))
      WinLose <- WinPost - LosePost
      WinDiff[num, 4:6] <- cbind(mean(WinLose), quantile(WinLose, .025), quantile(WinLose, .975))
    }
    num <- num + 1
  }
  return(data.frame(cbind(Seq, Est, High, Low, WinDiff)))
}

################################################################################
############ Figure 1: Male vs. Female Primary Elections by Party, 1972 - 2010
################################################################################

pdf(paste0(figure_dir, "MF_Races_ByParty.pdf"), height=7, width=10)
ggplot(dat, aes(factor(year), fill=factor(dem))) + 
  geom_bar(width=.75, position="dodge", alpha=.75) +
  xlab("Year") +
  ylab("Number of Male vs. Female House Primaries") + 
  scale_fill_manual(name="Party:", values = c("grey", "black"), labels=c("Republican", "Democrat")) +
  theme(legend.position="bottom") 
dev.off()

################################################################################
############ Figure 2: Female Candidate Primary Win Rates by Party, 1972 - 2010
################################################################################

###### Calculating Win Rates by Party
Agg_Rates <- data.frame(sort(unique(dat$year)), NA, NA, NA)
colnames(Agg_Rates) <- c("year", "femwin_perc", "femwin_perc_Rep", "femwin_perc_Dem")
for(i in unique(dat$year)){
  Agg_Rates[Agg_Rates$year == i,]$femwin_perc <- mean(dat[dat$year == i,]$primary_winner, na.rm=T)
  Agg_Rates[Agg_Rates$year == i,]$femwin_perc_Rep <- mean(dat[dat$year == i & dat$dem == 0,]$primary_winner, na.rm=T)
  Agg_Rates[Agg_Rates$year == i,]$femwin_perc_Dem <- mean(dat[dat$year == i & dat$dem == 1,]$primary_winner, na.rm=T)
}

#### Plot Rates by Party
pdf(paste0(figure_dir, "F_WinRates_ByParty.pdf"), height = 8, width = 12)
plot(x=Agg_Rates$year, y=Agg_Rates$femwin_perc_Rep, ylim=c(0, 1), xaxt = "n", pch=16, cex=.75, las=1, col = "darkred", 
     #main = "Primary Win Rates for Women Candidates by Party", 
     ylab="Primary Win Rate", xlab="Year")
axis(side = 1, at = seq(1972, 2010, 2))
lines(x=Agg_Rates$year, y=Agg_Rates$femwin_perc_Rep, ylim=c(0, 1), pch=16, cex=.75, col = "red", lty =2)
abline(h=.5, lty=3, col="gray")
## Dem
points(x=Agg_Rates$year, y=Agg_Rates$femwin_perc_Dem, ylim=c(0, 1), pch=16, cex=.75, las=1, col = "darkblue")
lines(x=Agg_Rates$year, y=Agg_Rates$femwin_perc_Dem, ylim=c(0, 1), pch=16, cex=.75, col = "blue")
legend("bottomright", c("Republican Women", "Democratic Women"), lty = c(2, 1), col = c("red", "blue"), cex = .75)
dev.off()

rm(i, Agg_Rates)


################################################################################
############ Table 1: Summary Statistics for Winning Candidates in `Male vs. Female' Primaries
################################################################################

## Adjusting Inc variable -- All rows are female candidates, so inc variable != variable of interest when male wins primary
dat$inc_for_summary <- ifelse(dat$primary_winner == 1, dat$inc, dat$Primary_Oppo_Inc)

## Calculating Group Means (Party by Gender)
out <- summaryBy(GenVote + GenWin + IncumbentOppo + candnumber + inc_for_summary + prez ~ party + primary_winner, data = dat)

## Creating/Exporting Table
out <- t(out)[2:8, ]
out <- apply(out, 2, function(x) round(as.numeric(x), 3))
out[1,] <- rep(c("Male Candidate", "Female Candidate"), 2)
colnames(out) <- c("Rep. Primary", "Rep. Primary", "Dem. Primary", "Dem. Primary")
row.names(out) <- c("Primary Winner", "General Election Vote Share", "General Election Win Percent",
                    "Nominated to Face Incumbent in General", "Number of Primary Challengers",
                    "Winning Candidate is Incumbent", "Presidential Vote Share, T-1")
addtorow <- list() ; addtorow$pos <- list(0, 7)
addtorow$command <- c(paste0(paste0('& \\multicolumn{2}{c}{', c("Republican Primaries", "Democratic Primaries"), '}', collapse=''), '\\\\'),
                      "\\hline & N = 257 & N = 227 & N = 351 & N = 407 \\\\")
table(dat$party, dat$primary_winner)
print(out)
print(xtable(out, align = "r|cc|cc", label = "Tab1", caption = "Summary Statistics for Winning Candidates in `Male vs. Female' Primaries"), 
      add.to.row=addtorow, 
      hline.after = c(1, 2, 3, 3, 4, 5, 6, 7, 7), 
      include.colnames = FALSE, 
      table.placement = "p!")

rm(out, addtorow)


################################################################################
############ Figure 3: General Election Vote Share Results
################################################################################

###### FIGURE 3 --- Saving as 2 x 2
pdf(paste0(figure_dir, "QuadPlot_VS.pdf"), height=10, width=10)
par(mfrow=c(2, 2))

###### Bandwidth Sequence:
Seq <- seq(.05, .25, .01)

#########
### Republicans
#########
Rep_Global <- lm(GenVote ~ primary_winner + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner, 
                 data=dat[dat$party == "R",])
summary(Rep_Global)

# Plot (Zoomed in)
BinPlot(Data = dat[dat$party == "R",], 
        RegObj = Rep_Global, 
        binwidth = .025, bandwidth = 1,
        runningVar = dat[dat$party == "R",]$PrimMargin,
        OutcomeCol = which(colnames(dat) == "GenVote"), 
        xlab="Primary Margin of Victory/Loss", ylab="General Election Vote Share", 
        title="Republican Primaries: Quadratic Polynomial", 
        color="black", xlim=c(-.3, .3), ylim=c(.2,.8))  
Predict.Plot(Rep_Global, pred.var = "PrimMargin", primary_winner = 0, add=T, data=dat[dat$party == "R" & dat$PrimMargin <= 0,], plot.args = list(lty=2, lwd=1.5, col='black'))
Predict.Plot(Rep_Global, pred.var = "PrimMargin", primary_winner = 1, add=T, data=dat[dat$party == "R" & dat$PrimMargin >= 0,], plot.args = list(lty=2, lwd=1.5, col='black'))

##Local Linear
Form <- formula(GenVote ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + inc + Y82to90 + Y92to00 + Y02to10)
RepVote <- LocalLin(Seq=Seq, Sub=dat$party == "R", Formula=Form, PredProb=F)
plotCI(x=RepVote$Seq, y=RepVote$Est, ui= RepVote$High, li= RepVote$Low, pch=16, cex=.8, 
       col=alpha("black", .75), las=1, ylim=c(-.3,.3),
       main="Republican Primaries: Local Linear", xlab="Bandwidth", 
       ylab=expression(beta ~ 1)); abline(h=0, lty=2)

#########
### Democrats
#########
Dem_Global <- lm(GenVote ~ primary_winner + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner, 
                 data=dat[dat$party == "D",])
summary(Dem_Global)

# Plot (Zoomed In)
BinPlot(Data = dat[dat$party %in% "D",], 
        RegObj = Dem_Global, 
        binwidth = .025, bandwidth = 1,
        runningVar = dat[dat$party %in% "D",]$PrimMargin,
        OutcomeCol = which(colnames(dat) == "GenVote"), 
        xlab="Primary Margin of Victory/Loss", ylab="General Election Vote Share", 
        title="Democratic Primaries: Quadratic Polynomial", 
        color="black", xlim=c(-.3, .3), ylim=c(.2,.8))  
Predict.Plot(Dem_Global, pred.var = "PrimMargin", primary_winner = 0, add=T, data=dat[dat$party == "D" & dat$PrimMargin <= 0,], plot.args = list(lty=2, lwd=1.5, col='black'))
Predict.Plot(Dem_Global, pred.var = "PrimMargin", primary_winner = 1, add=T, data=dat[dat$party == "D" & dat$PrimMargin >= 0,], plot.args = list(lty=2, lwd=1.5, col='black'))

### Local Linear
Form <- formula(GenVote ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + inc + Y82to90 + Y92to00 + Y02to10)
DemVote <- LocalLin(Seq=Seq, Sub=dat$party == "D", Formula=Form, PredProb=F)
plotCI(x=DemVote$Seq, y=DemVote$Est, ui= DemVote$High, li= DemVote$Low, pch=16, cex=.8, 
       col=alpha("black", .75), las=1,ylim=c(-.3,.3), main="Democratic Primaries: Local Linear", 
       xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)


dev.off()
### ** End Plot PDF ** 


################################################################################
############ Tables 3 & 4: General Election Vote Share Results
################################################################################

Form <- formula(GenVote ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + inc + Y82to90 + Y92to00 + Y02to10)
BW <- .2  ; BW.5 <- .1  ;  BW.25 <- .05     # Bandwidths
Rep_LL <- lm(Form, data=dat[dat$party == "R" & abs(dat$PrimMargin) <= BW/2,])
Rep_LL.5 <- lm(Form, data=dat[dat$party == "R" & abs(dat$PrimMargin) <= BW.5/2,])
Rep_LL.25 <- lm(Form, data=dat[dat$party == "R" & abs(dat$PrimMargin) <= BW.25/2,])
Dem_LL <- lm(Form, data=dat[dat$party == "D" & abs(dat$PrimMargin) <= BW/2,])
Dem_LL.5 <- lm(Form, data=dat[dat$party == "D" & abs(dat$PrimMargin) <= BW.5/2,])
Dem_LL.25 <- lm(Form, data=dat[dat$party == "D" & abs(dat$PrimMargin) <= BW.25/2,])

### Using Max of Conventional/Robust SEs
SE_Func <- function(Model){
  Est <- coef(Model)
  RegSE <- summary(Model)$coefficients[, 2] #Standard SEs
  Robust <- coeftest(Model, vcov = hccm(Model, type="hc1"))[,2] #Robust SE's
  OutSE <- apply(data.frame(RegSE, Robust), 1, function(x) max(x))
  return(OutSE)
}

Rep_SE <- SE_Func(Rep_LL) 
Rep_SE.5 <- SE_Func(Rep_LL.5) 
Rep_SE.25 <- SE_Func(Rep_LL.25) 
Dem_SE <- SE_Func(Dem_LL) 
Dem_SE.5 <- SE_Func(Dem_LL.5) 
Dem_SE.25 <- SE_Func(Dem_LL.25) 

# With Covariates
Rep_Globalc <- lm(GenVote ~ primary_winner + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner + prez + PrezYear + inc + Y82to90 + Y92to00 + Y02to10, data=dat[dat$party == "R",])
Dem_Globalc <- lm(GenVote ~ primary_winner + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner + prez + PrezYear + inc + Y82to90 + Y92to00 + Y02to10, data=dat[dat$party == "D",])

# Table 3: Polynomial VS Models with Covariates
# stargazer(Rep_Global, Rep_Globalc, Dem_Global, Dem_Globalc, omit.stat=c("adj.rsq", "f", "ser"), type = "text")
stargazer(Rep_Global, Rep_Globalc, Dem_Global, Dem_Globalc,
          omit.stat=c("adj.rsq", "f", "ser"),
          dep.var.labels="General Election Vote Share", 
          #label="VS_poly", 
          style="apsr", 
          star.char = c("+", "*", "**"),
          title="Quadratic Polynomial Regressions of General Election Vote Share on Female Candidate Primary Win",
          covariate.labels=c("FemaleWin", "Margin", "Margin$^2$", "Pres. Vote T-1", "Pres. Elec. Year", "Incumbent", 
                             "1982 - 1990", "1992 - 2000", "2002 - 2010", "FemaleWin*Margin", "FemaleWin*Margin$^2$"),
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01", notes.append = FALSE)

#### Table 4: Local Linear VS Models
# stargazer(Rep_LL, Rep_LL.5, Rep_LL.25, Dem_LL, Dem_LL.5, Dem_LL.25, omit.stat=c("adj.rsq", "f", "ser"), type = "text")
stargazer(Rep_LL, Rep_LL.5, Rep_LL.25, Dem_LL, Dem_LL.5, Dem_LL.25, 
          omit.stat=c("adj.rsq", "f", "ser"), 
          dep.var.labels="General Election Vote Share",
          #label="VS_LL", 
          style="apsr", 
          star.char = c("+", "*", "**"),
          title="Local Linear Regressions of General Election Vote Share on Female Candidate Primary Win",
          covariate.labels=c("FemaleWin", "Margin", "Pres. Vote T-1", "Pres. Elec. Year", "Incumbent", 
                             "1982 - 1990", "1992 - 2000", "2002 - 2010", "FemaleWin*Margin"),
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01", notes.append = FALSE,
          se=list(Rep_SE, Rep_SE.5, Rep_SE.25, Dem_SE, Dem_SE.5, Dem_SE.25), p.auto=T)

rm(Rep_Global, Rep_Globalc, Dem_Global, Dem_Globalc, BW, BW.25, BW.5, Seq, DemVote, RepVote, Form,
   Rep_LL, Rep_LL.5, Rep_LL.25, Dem_LL, Dem_LL.5, Dem_LL.25, Rep_SE, Rep_SE.5, Rep_SE.25, Dem_SE, Dem_SE.5, Dem_SE.25)


################################################################################
############ Figure 4: Win Probability Results
################################################################################

###### Bandwidth Sequence:
Seq <- seq(.05, .25, .01)

###### Saving as 2 x 2 Figure
pdf(paste0(figure_dir, "QuadPlot_Win.pdf"), height=10, width=10)
par(mfrow=c(2, 2))

#########
### Republicans
#########

## Polynomial Win Probability Models (With and Without Covariates)
Rep_Global_W <- lm(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner, data=dat[dat$party == "R",])
summary(Rep_Global_W)

Rep_Global_W2 <- lm(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner + prez + PrezYear + inc + Y82to90 + Y92to00 + Y02to10, data=dat[dat$party == "R",])
summary(Rep_Global_W2)

### Note: Findings do not change with Year Fixed Effects, with or without Covariates
summary(lm(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner + factor(year), data=dat[dat$party == "R",]))
summary(lm(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner + prez + PrezYear + inc + factor(year), data=dat[dat$party == "R",]))

### Plot 
BinPlot(Data = dat[dat$party == "R",], 
        RegObj = Rep_Global_W, 
        binwidth = .025, bandwidth = 1,
        runningVar = dat[dat$party == "R",]$PrimMargin,
        OutcomeCol = which(colnames(dat) == "GenWin"), 
        xlab="Primary Margin of Victory/Loss", ylab="General Election Win", 
        title="Republican Primaries: Quadratic Polynomial", 
        color="black", xlim=c(-.3,.3), ylim=c(0,1))     
Predict.Plot(Rep_Global_W, pred.var = "PrimMargin", primary_winner = 0, add=T, data=dat[dat$party == "R" & dat$PrimMargin <= 0,], plot.args = list(lty=2, lwd=1.5, col='black'))
Predict.Plot(Rep_Global_W, pred.var = "PrimMargin", primary_winner = 1, add=T, data=dat[dat$party == "R" & dat$PrimMargin >= 0,], plot.args = list(lty=2, lwd=1.5, col='black'))

##Local Linear
Form <- formula(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + inc + Y82to90 + Y92to00 + Y02to10)
RepWin <- LocalLin(Seq=Seq, Sub=dat$party == "R", Formula=Form, PredProb=F)
plotCI(x=RepWin$Seq, y=RepWin$Est, ui= RepWin$High, li= RepWin$Low, pch=16, cex=.8, col=alpha("black", .75), las=1, ylim=c(-1, .35),
       main="Republican Primaries: Local Linear", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

##########
### Democrats
##########

Dem_Global_W <- lm(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner, data=dat[dat$party == "D",])
Dem_Global_W2 <- lm(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner + prez + PrezYear + inc + Y82to90 + Y92to00 + Y02to10, data=dat[dat$party == "D",])
summary(Dem_Global_W)

BinPlot(Data = dat[dat$party == "D",], 
        RegObj = Dem_Global_W, 
        binwidth = .025, bandwidth = 1,
        runningVar = dat[dat$party == "D",]$PrimMargin,
        OutcomeCol = which(colnames(dat) == "GenWin"), 
        xlab="Primary Margin of Victory/Loss", ylab="General Election Win", 
        title="Democratic Primaries: Quadratic Polynomial", 
        color="black", xlim=c(-.3,.3), ylim=c(0,1))  
Predict.Plot(Dem_Global_W, pred.var = "PrimMargin", primary_winner = 0, add=T, data=dat[dat$party == "D" & dat$PrimMargin <= 0,], plot.args = list(lty=2, lwd=1.5, col='black'))
Predict.Plot(Dem_Global_W, pred.var = "PrimMargin", primary_winner = 1, add=T, data=dat[dat$party == "D" & dat$PrimMargin >= 0,], plot.args = list(lty=2, lwd=1.5, col='black'))

### Local Linear
DemWin <- LocalLin(Seq=Seq, Sub=dat$party == "D", Formula=Form, PredProb=F)
plotCI(x=DemWin$Seq, y=DemWin$Est, ui= DemWin$High, li= DemWin$Low, pch=16, cex=.8, col=alpha("black", .75),las=1, ylim=c(-1, .35),
       main="Democratic Primaries: Local Linear", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

dev.off()
### *** END PLOT ***

########################
#### Footnote 18 - Testing Null of Equality of Coefficients
########################

### Quadratic Model
Form <- formula(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + primary_winner*I(PrimMargin^2))
R_test <- lm(Form, data=dat[dat$party == "R",])
D_test <- lm(Form, data=dat[dat$party == "D",])
R_out <- rmvnorm(10000, mean=coef(R_test), sigma=vcov(R_test))
D_out <- rmvnorm(10000, mean=coef(D_test), sigma=vcov(D_test))
quantile(R_out[,2] - D_out[,2], c(.025, .975))

### Local Linear Model
Form <- formula(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + inc + Y82to90 + Y92to00 + Y02to10)
for(i in c(10:30)){
  bw <- i/100
  R_test <- lm(Form, data=dat[dat$party == "R" & abs(dat$PrimMargin) <= bw/2,])
  D_test <- lm(Form, data=dat[dat$party == "D" & abs(dat$PrimMargin) <= bw/2,])
  R_out <- rmvnorm(10000, mean=coef(R_test), sigma=vcov(R_test))
  D_out <- rmvnorm(10000, mean=coef(D_test), sigma=vcov(D_test))
  Stat <- R_out[,2] - D_out[,2]
  Q <- round(quantile(Stat, c(.025, .975)), 3)
  print(paste("Bandwidth: ", bw, " --- Confidence Interval: ", Q[1], "-", Q[2], sep=""))
}

rm(Form, bw, R_test, D_test, R_out, D_out, Stat, Q, i)

################################################################
### Table 5 + 6: RDD Results - General Election Win on Female Candidate Primary Win
################################################################

Form <- formula(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + inc + Y82to90 + Y92to00 + Y02to10)
BW <- .20  ;  BW.5 <- .10  ;  BW.25 <- .05 

Rep_LL <- lm(Form, data=dat[dat$party == "R" & abs(dat$PrimMargin) <= BW/2,])
Rep_LL.5 <- lm(Form, data=dat[dat$party == "R" & abs(dat$PrimMargin) <= BW.5/2,])
Rep_LL.25 <- lm(Form, data=dat[dat$party == "R" & abs(dat$PrimMargin) <= BW.25/2,])
Dem_LL <- lm(Form, data=dat[dat$party == "D" & abs(dat$PrimMargin) <= BW/2,])
Dem_LL.5 <- lm(Form, data=dat[dat$party == "D" & abs(dat$PrimMargin) <= BW.5/2,])
Dem_LL.25 <- lm(Form, data=dat[dat$party == "D" & abs(dat$PrimMargin) <= BW.25/2,])

### Max of Conventional and Robust SE's
Rep_SE <- SE_Func(Rep_LL) 
Rep_SE.5 <- SE_Func(Rep_LL.5) 
Rep_SE.25 <- SE_Func(Rep_LL.25)  
Dem_SE <- SE_Func(Dem_LL) 
Dem_SE.5 <- SE_Func(Dem_LL.5)  
Dem_SE.25 <- SE_Func(Dem_LL.25) 

#### Table 5: Polynomial Win Models
# stargazer(Rep_Global_W, Rep_Global_W2, Dem_Global_W, Dem_Global_W2, omit.stat=c("adj.rsq", "f", "ser"), type = "text")
stargazer(Rep_Global_W, Rep_Global_W2, Dem_Global_W, Dem_Global_W2, 
          omit.stat=c("adj.rsq", "f", "ser"), 
          dep.var.labels="General Election Win", 
          style="apsr", 
          #label="Win_poly", 
          star.char = c("+", "*", "**"), 
          title="Quadratic Polynomial Regressions of General Election Win on Female Candidate Primary Win",
          covariate.labels=c("FemaleWin", "Margin", "Margin$^2$", "Pres. Vote T-1", "Pres. Elec. Year", "Incumbent", 
                             "1982 - 1990", "1992 - 2000", "2002 - 2010", "FemaleWin*Margin", "FemaleWin*Margin$^2$"),
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01", notes.append = FALSE)

#### Table 6: Local Linear Win Models
# stargazer(Rep_LL, Rep_LL.5, Rep_LL.25, Dem_LL, Dem_LL.5, Dem_LL.25, omit.stat=c("adj.rsq", "f", "ser"), type = "text")
stargazer(Rep_LL, Rep_LL.5, Rep_LL.25, Dem_LL, Dem_LL.5, Dem_LL.25, 
          omit.stat=c("adj.rsq", "f", "ser"), 
          dep.var.labels="General Election Win", 
          style="apsr", 
          #label="Win_LL", 
          star.char = c("+", "*", "**"), 
          title="Local Linear Regressions of General Election Win on Female Candidate Primary Win",
          covariate.labels=c("FemaleWin", "Margin", "Pres. Vote T-1", "Pres. Elec. Year", "Incumbent", 
                             "1982 - 1990", "1992 - 2000", "2002 - 2010", "FemaleWin*Margin"),
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01", notes.append = FALSE,
          se=list(Rep_SE, Rep_SE.5, Rep_SE.25, Dem_SE, Dem_SE.5, Dem_SE.25), p.auto=T)

rm(Rep_Global_W, Rep_Global_W2, Dem_Global_W, Dem_Global_W2, BW, BW.5, BW.25, SE_Func, RepWin, DemWin,
   Rep_LL, Rep_LL.5, Rep_LL.25, Dem_LL, Dem_LL.5, Dem_LL.25, Rep_SE, Rep_SE.5, Rep_SE.25, Dem_SE, Dem_SE.5, Dem_SE.25)

################################################################################
############ Figure 5: Variance in Outcomes - Male versus Female
################################################################################

dat$Win_Factor <- factor(dat$primary_winner)

## Full Distribution
full <- ggplot(dat[dat$party == "R",], aes(x=GenVote, ..density.., col = Win_Factor, linetype = Win_Factor)) + 
  geom_freqpoly(position='identity', binwidth = .05) + 
  scale_linetype_manual(values = c("solid", "dashed"), 
                        guide = guide_legend(title=""),
                        labels = c("Male Candidate Wins", "Female Candidate Wins")) +
  scale_color_manual(values = c("grey50", "black"), 
                     guide = guide_legend(title=""),
                     labels = c("Male Candidate Wins", "Female Candidate Wins")) + 
  theme_bw() + 
  theme(legend.position="bottom", 
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line = element_line(colour = "black")) + 
  xlab("General Election Vote Share") + 
  ylab("Density") + ylim(0, 5) +
  ggtitle("Full Sample") + 
  geom_vline(xintercept = mean(dat[dat$party == "R" & dat$primary_winner == 0,]$GenVote), lty = 1, lwd = .35, col = "grey50") + 
  geom_vline(xintercept = mean(dat[dat$party == "R" & dat$primary_winner == 1,]$GenVote), lty = 2, lwd = .35, col = "black")  

## Distribution for Close Races = 15% Bandwidth (+/- 7.5%)
sub <- ggplot(dat[dat$party == "R" & abs(dat$PrimMargin) < .075,], aes(x=GenVote, ..density.., col = Win_Factor, linetype = Win_Factor)) + 
  geom_freqpoly(position='identity', binwidth = .05) + # for density: aes(y=..density..),
  scale_linetype_manual(values = c("solid", "dashed"), guide = guide_legend(title=""),
                        labels = c("Male Candidate Wins", "Female Candidate Wins")) +
  scale_color_manual(values = c("grey50", "black"), guide = guide_legend(title=""),
                     labels = c("Male Candidate Wins", "Female Candidate Wins")) + 
  theme_bw() + 
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(), 
        axis.line = element_line(colour = "black")) + 
  xlab("General Election Vote Share") + 
  ylab("Density") + ylim(0, 5) + 
  ggtitle("15 Percent Bandwidth") + 
  geom_vline(xintercept = mean(dat[dat$party == "R" & dat$primary_winner == 0,]$GenVote), lty = 1, lwd = .35, col = "grey50") + 
  geom_vline(xintercept = mean(dat[dat$party == "R" & dat$primary_winner == 1,]$GenVote), lty = 2, lwd = .35, col = "black")

## Extract legend -- via http://stackoverflow.com/questions/13649473/add-a-common-legend-for-combined-ggplots
g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}
mylegend <- g_legend(full)

## Combined Plots + Save
pdf(paste0(figure_dir, "Rep_GenVote_Dist_Comparison.pdf"), height = 5, width = 10)
grid.arrange(arrangeGrob(full + theme(legend.position="none"), 
                         sub + theme(legend.position="none"), nrow=1),
             mylegend, 
             nrow=2, 
             #top = "Distribution of General Election Vote Share for Republican Candidates",
             heights=c(6, 1))
dev.off()

rm(full, sub, mylegend, g_legend)

####################
#### F-Test for Ratio of Variances (p. 17, Ungated Version)
####################

### 5% Bandwidth (+/- 2.5%)
var.test(dat[dat$party == "R" & abs(dat$PrimMargin) < .025 & dat$primary_winner == 0,]$GenVote, 
         dat[dat$party == "R" & abs(dat$PrimMargin) < .025  & dat$primary_winner == 1,]$GenVote)

### 10% Bandwidth (+/- 5%)
var.test(dat[dat$party == "R" & abs(dat$PrimMargin) < .05 & dat$primary_winner == 0,]$GenVote, 
         dat[dat$party == "R" & abs(dat$PrimMargin) < .05  & dat$primary_winner == 1,]$GenVote)

### 15% Bandwidth (+/- 7.5%)
var.test(dat[dat$party == "R" & abs(dat$PrimMargin) < .075 & dat$primary_winner == 0,]$GenVote, 
         dat[dat$party == "R" & abs(dat$PrimMargin) < .075  & dat$primary_winner == 1,]$GenVote)

### Full Sample
var.test(dat[dat$party == "R" & dat$primary_winner == 0,]$GenVote, 
         dat[dat$party == "R" & dat$primary_winner == 1,]$GenVote)


################################################################################
############ Figure 6: Pre/Post 1990 - Vote Share
################################################################################

### Adjusting Formula --- Insufficient Inc Variation within interaction groups
Form <- formula(GenVote ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + inc + primary_winner*Post90)
Form_no_inc <- formula(GenVote ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + primary_winner*Post90)
Seq <- seq(.05, .25, .01) 

########
## Republicans 
Rep_VS_Diff <- LocalLin(Seq=Seq, Sub=dat$party == "R", Formula=Form, PredProb=T, use_alt = T, alt_Form = Form_no_inc, alt_bw = .09)

#######
## Democrats 
Dem_VS_Diff <- LocalLin(Seq=Seq, Sub=dat$party == "D", Formula=Form, PredProb=T, use_alt = F, alt_Form = NA, alt_bw = NA)

################
## Figure 6 (See LocalLin Function at start for details on estimates/CIs):
pdf(paste0(figure_dir, "VS_PrePost.pdf"), height=10, width=10)
par(mfrow=c(2,2))

#### Republicans
# Predicted Difference in Vote Share for Female Republican Candidate Nominated After - Before 1990
plotCI(Rep_VS_Diff$Seq, y=Rep_VS_Diff$V5, ui= Rep_VS_Diff$V7, li= Rep_VS_Diff$V6, pch=16, cex=.8, ylim=c(-.25, .25), 
       col=alpha("black", .75), main="Female Republicans, Before and After 1990", slty=5,
       xlab="Bandwidth", ylab="Predicted Difference in Vote Share"); abline(h=0, lty=2)
# Predicted Difference in Vote Share for Female Republican Candidate (Post 1990) - Male Republican Candidate (Post 1990)
plotCI(Rep_VS_Diff$Seq, y=Rep_VS_Diff$V8, ui= Rep_VS_Diff$V10, li= Rep_VS_Diff$V9, pch=16, cex=.8, ylim=c(-.25, .25),
       col=alpha("black", .75), main="Female vs. Male Republicans, After 1990", slty=5,
       xlab="Bandwidth", ylab="Predicted Difference in Vote Share"); abline(h=0, lty=2)

#### Democrats
#Predicted Difference in Vote Share for Female Democratic Candidate Nominated After - Before 1990
plotCI(Dem_VS_Diff$Seq, y=Dem_VS_Diff$V5, ui= Dem_VS_Diff$V7, li= Dem_VS_Diff$V6, pch=16, cex=.8, ylim=c(-.25, .25),
       col=alpha("black", .75), main="Female Democrats, Before and After 1990", 
       xlab="Bandwidth", ylab="Predicted Difference in Vote Share"); abline(h=0, lty=2)
#Predicted Difference in Vote Share for Female Democratic Candidate (Post 1990) - Male Democratic Candidate (Post 1990)
plotCI(Dem_VS_Diff$Seq, y=Dem_VS_Diff$V8, ui= Dem_VS_Diff$V10, li= Dem_VS_Diff$V9, pch=16, cex=.8,  ylim=c(-.25, .25),
       col=alpha("black", .75), main="Female vs. Male Democrats, After 1990",
       xlab="Bandwidth", ylab="Predicted Difference in Vote Share"); abline(h=0, lty=2)
dev.off()
## ** END PLOT **


################################################################################
############ Figure 7: Pre/Post 1990 - Win Probability
################################################################################

### Adjusting Formula --- Insufficient Inc Variation within interaction groups
Form <- formula(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + inc + primary_winner*Post90)
Form_no_inc <- formula(GenWin ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + primary_winner*Post90)
Seq <- seq(.05, .25, .01) 

########
## Republicans 
Rep_Win_Diff <- LocalLin(Seq=Seq, Sub=dat$party == "R", Formula=Form, PredProb=T, use_alt = T, alt_Form = Form_no_inc, alt_bw = .09)

#######
## Democrats 
Dem_Win_Diff <- LocalLin(Seq=Seq, Sub=dat$party == "D", Formula=Form, PredProb=T, use_alt = F)

################
## Figure 7 
pdf(paste0(figure_dir, "Win_PrePost.pdf"), height=10, width=10)
par(mfrow=c(2,2))

#### Republicans
# Predicted Difference in Win Probability for Female Republican Candidate Nominated After - Before 1990
plotCI(Rep_Win_Diff$Seq, y=Rep_Win_Diff$V5, ui= Rep_Win_Diff$V7, li= Rep_Win_Diff$V6, pch=16, cex=.8, ylim=c(-1, .5),
       col=alpha("black", .75), main="Female Republicans: Before and After 1990", slty=5,
       xlab="Bandwidth", ylab="Predicted Difference in Win Probability"); abline(h=0, lty=2)
# Predicted Difference in Win Prob. for Female Republican Candidate (Post 1990) - Male Republican Candidate (Post 1990)
plotCI(Rep_Win_Diff$Seq, y=Rep_Win_Diff$V8, ui= Rep_Win_Diff$V10, li= Rep_Win_Diff$V9, pch=16, cex=.8,  ylim=c(-1, .5),
       col=alpha("black", .75), main="Female vs. Male Republicans, After 1990", slty=5,
       xlab="Bandwidth", ylab="Predicted Difference in Win Probability"); abline(h=0, lty=2)

#### Democrats
# Predicted Difference in Win Probability for Female Democratic Candidate Nominated After - Before 1990
plotCI(Dem_Win_Diff$Seq, y=Dem_Win_Diff$V5, ui= Dem_Win_Diff$V7, li= Dem_Win_Diff$V6, pch=16, cex=.8,  ylim=c(-1, .5),
       col=alpha("black", .75), main="Female Democrats, Before and After 1990",
       xlab="Bandwidth", ylab="Predicted Difference in Win Probability"); abline(h=0, lty=2)
# Predicted Difference in Win Prob. for Female Democratic Candidate (Post 1990) - Male Democratic Candidate (Post 1990)
plotCI(Dem_Win_Diff$Seq, y=Dem_Win_Diff$V8, ui= Dem_Win_Diff$V10, li= Dem_Win_Diff$V9, pch=16, cex=.8,  ylim=c(-1, .5),
       col=alpha("black", .75), main="Female vs. Male Democrats, After 1990",
       xlab="Bandwidth", ylab="Predicted Difference in Win Probability"); abline(h=0, lty=2)

dev.off()
## ** END PLOT **

rm(Dem_VS_Diff, Dem_Win_Diff, Rep_VS_Diff, Rep_Win_Diff, Form, Seq, Form_no_inc)

############################################
### Figure 9 (APPENDIX) - Year Interaction Specification
#############################################

time_inter <- lm(GenWin ~ primary_winner*year + PrimMargin + primary_winner*PrimMargin + primary_winner*I(PrimMargin^2), 
                 data = dat[dat$party == "R",])
time_inter2 <- lm(GenWin ~ primary_winner*year + PrimMargin + primary_winner*PrimMargin + primary_winner*I(PrimMargin^2), 
                 data = dat[dat$party == "R" & abs(dat$PrimMargin) < .20,])
time_inter3 <- lm(GenWin ~ primary_winner*year + PrimMargin + primary_winner*PrimMargin + primary_winner*I(PrimMargin^2) + 
                    prez + PrezYear + inc , data = dat[dat$party == "R",])

# Estimate of coefficient on primary_winner over time
time1 <- interplot(time_inter, "primary_winner", "year", hist = TRUE) + geom_hline(aes(yintercept = 0), lty = 2, lwd = .5)
time2 <- interplot(time_inter2, "primary_winner", "year", hist = TRUE) + geom_hline(aes(yintercept = 0), lty = 2, lwd = .5)
time3 <- interplot(time_inter3, "primary_winner", "year", hist = TRUE) + geom_hline(aes(yintercept = 0), lty = 2, lwd = .5)

time1 <- time1 + xlab("Year") + ylab("General Election Win Probability") + ggtitle("Base Model (No Controls)")
time2 <- time2 + xlab("Year") + ylab("General Election Win Probability") + ggtitle("20 Percent Bandwidth (No Controls)")
time3 <- time3 + xlab("Year") + ylab("General Election Win Probability") + ggtitle("Full Model (With Controls)")

pdf(paste0(figure_dir, "Win_Time_Interact.pdf"), heigh=4, width=12)
grid.arrange(time1, time2, time3, 
             #top = "Variation in Effect on General Election Win Probability across Election Cycles (Republicans)"
             nrow = 1)
dev.off()

rm(time1, time2, time3, time_inter, time_inter2, time_inter3)


################################################################################
############ Campaign Finance Results
################################################################################

### Subsetting to relevant data
cf_data <- dat[dat$year >= 1990,]

#Decade Indicators
cf_data$Y90s <- ifelse(cf_data$year >= 1990 & cf_data$year <= 1998, 1, 0)
cf_data$Y2000s <- ifelse(cf_data$year >= 2000 & cf_data$year <= 2010, 1, 0)

## Post Bipartisan Campaign Reform Act of 2002
cf_data$BCRA <- ifelse(cf_data$year > 2002, 1, 0)

### Share of General Election (and Primary) Contributions Received
cf_data$GenShare <- cf_data$WinFunds/(cf_data$GenOppoGenFunds + cf_data$WinFunds) * 100
cf_data$PrimShare <- cf_data$PrimaryFunds/(cf_data$OppoPrimFunds + cf_data$PrimaryFunds) * 100

## Adjusting Primary Contributions to 0 if no record
cf_data$OppoPrimFunds2 <- ifelse(cf_data$OppoPrimFunds %in% NA, 0, cf_data$OppoPrimFunds)
cf_data$PrimaryFunds2 <- ifelse(cf_data$PrimaryFunds %in% NA, 0, cf_data$PrimaryFunds)
cf_data$PrimShare_adj <- ifelse(cf_data$PrimaryFunds2 == 0 & cf_data$OppoPrimFunds2 == 0, NA, 100 * cf_data$PrimaryFunds2 / (cf_data$PrimaryFunds2 + cf_data$OppoPrimFunds2))

################################################################################
############ Table 2: Regressions with Campaign Finance DV
################################################################################

####### Republicans 
RepShare <- lm(GenShare ~ primary_winner + PrimMargin + I(primary_winner * PrimMargin)  + I(PrimMargin^2) + I(primary_winner * PrimMargin^2), 
               data=cf_data[cf_data$party == "R",])
RepShare2 <- lm(GenShare ~ primary_winner + PrimMargin + I(PrimMargin^2) + I(primary_winner * PrimMargin) + I(primary_winner * PrimMargin^2) + 
                  prez + inc + BCRA, data=cf_data[cf_data$party == "R",])
RepShare3 <- lm(GenShare ~ primary_winner + PrimMargin + I(PrimMargin^2) + I(primary_winner * PrimMargin) + I(primary_winner * PrimMargin^2) + 
                  prez + inc + BCRA + PrimShare_adj, data=cf_data[cf_data$party == "R",])

##### Democrats
DemShare <- lm(GenShare ~ primary_winner + PrimMargin + I(primary_winner * PrimMargin)  + I(PrimMargin^2) + I(primary_winner * PrimMargin^2), 
               data=cf_data[cf_data$party == "D",])
DemShare2 <- lm(GenShare ~ primary_winner + PrimMargin + I(PrimMargin^2) + I(primary_winner * PrimMargin) + I(primary_winner * PrimMargin^2) + 
                  prez + inc + BCRA, data=cf_data[cf_data$party == "D",])
DemShare3 <- lm(GenShare ~ primary_winner + PrimMargin + I(PrimMargin^2) + I(primary_winner * PrimMargin) + I(primary_winner * PrimMargin^2) + 
                  prez + inc + BCRA + PrimShare_adj, data=cf_data[cf_data$party == "D",])

#### Table 2:
# stargazer(RepShare, RepShare2, RepShare3, DemShare, DemShare2, DemShare3, omit.stat=c("adj.rsq", "f", "ser"), type = "text")
stargazer(RepShare, RepShare2, RepShare3, DemShare, DemShare2, DemShare3,
          dep.var.labels="General Election Contribution Share", 
          label="share",
          omit.stat=c("adj.rsq", "f", "ser"), 
          omit="year", 
          star.char = c("+", "*", "**"),
          covariate.labels=c("FemaleWin", "Margin", "FemaleWin * Margin", "Margin$^2$", "FemaleWin*Margin$^2$", 
                             "Pres. Vote T-1", "Incumbent", "Post-BCRA", "Primary Share", "Constant"),
          notes = "$^{+}$p$<$0.1; $^{*}$p$<$0.05; $^{**}$p$<$0.01", notes.append = FALSE,
          title="Quadratic Polynomial Regressions of General Election Contribution Share on Female Candidate Primary Win")

rm(RepShare, RepShare2, RepShare3, DemShare, DemShare2, DemShare3)

###########################################################
############ Figure 8: Campaign Finance Interaction - Time * Female Candidate Priamry Win
############################################################

# Republicans
rep_interact_cash <- lm(GenShare ~ primary_winner*year + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner, 
                        data=cf_data[cf_data$party == "R",])
rep_interact_cash2 <- lm(GenShare ~ primary_winner*year + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner + prez + inc + BCRA, 
                         data=cf_data[cf_data$party == "R",])
rep_interact_cash3 <- lm(GenShare ~ primary_winner*year + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner + prez + inc + BCRA + PrimShare_adj, 
                         data=cf_data[cf_data$party == "R",])

time1 <- interplot(rep_interact_cash, "primary_winner", "year", hist = TRUE) + geom_hline(aes(yintercept = 0), lty = 2, lwd = .5)
time2 <- interplot(rep_interact_cash2, "primary_winner", "year", hist = TRUE) + geom_hline(aes(yintercept = 0), lty = 2, lwd = .5)
time3 <- interplot(rep_interact_cash3, "primary_winner", "year", hist = TRUE) + geom_hline(aes(yintercept = 0), lty = 2, lwd = .5)

time1 <- time1 + xlab("Year") + ylab("General Election Contribution Share") + ggtitle("Base Model")
time2 <- time2 + xlab("Year") + ylab("General Election Contribution Share") + ggtitle("Partial Controls")
time3 <- time3 + xlab("Year") + ylab("General Election Contribution Share") + ggtitle("Full Controls")

pdf(paste0(figure_dir, "CF_Time_Interact.pdf"), heigh=4, width=12)
grid.arrange(time1, time2, time3, 
             #top = "Variation in Effect on Contribution Share acrss Election Cycles (Republicans)",
             nrow = 1)
dev.off()

rm(time1, time2, time3, rep_interact_cash, rep_interact_cash2, rep_interact_cash3)

###########################
#### Figure 10 - Campaign Finance, Local Linear 
###########################

pdf(paste0(figure_dir, "CF_LocalLin.pdf"), heigh=5, width=10)
par(mfrow=c(1,2))

Coef <- rep(NA,21) ; SE <- rep(NA, 21) ; num <- 1
for(i in seq(.05, .25, by=.01)){
  Local <- lm(GenShare ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + inc + BCRA, data=cf_data, 
              subset=cf_data$party == "R" & abs(cf_data$PrimMargin) <= i)
  Coef[num] <- coef(Local)[2]
  Reg <- summary(Local)$coefficients[2, 2] #Standard SEs
  VarCovar <- hccm(Local, type="hc1")
  Robust <- coeftest(Local, vcov = VarCovar)[2,2] #Robust SE's
  SE[num] <- max(Reg, Robust) #Using maximum of regular and Robust SE's
  num <- num + 1
}
####### Including adjusted share variable has minimal effects on LL estimates

#Plotting the Variation in coefficients by  bandwidth size
Upper <- Coef + 1.96*SE ; Lower <- Coef - 1.96*SE

X <- seq(.05, .25, by=.01)
plotCI(x=X, y=Coef, ui= Upper, li= Lower, pch=16, cex=.8, ylim=c(-85, 20), ylab=expression(beta ~ 1), 
       xlab="Bandwidth Size", main="Republican Primaries") #Local Linear RDDs
abline(h=0, lty=2)

###############
### DEM

Coef <- rep(NA,21) ; SE <- rep(NA, 21) ; num <- 1
for(i in seq(.05, .25, by=.01)){
  Local <- lm(GenShare ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + inc + BCRA, 
              data=cf_data, subset=cf_data$party == "D" & abs(cf_data$PrimMargin) <= i)
  Coef[num] <- coef(Local)[2]
  Reg <- summary(Local)$coefficients[2, 2] #Standard SEs
  VarCovar <- hccm(Local, type="hc1")
  Robust <- coeftest(Local, vcov = VarCovar)[2,2] #Robust SE's
  SE[num] <- max(Reg, Robust) #Using maximum of regular and Robust SE's
  num <- num + 1
}

#Plotting the Variation in coefficients by  bandwidth size
Upper <- Coef + 1.96*SE ; Lower <- Coef - 1.96*SE; X <- seq(.05, .25, by=.01)
plotCI(x=X, y=Coef, ui= Upper, li= Lower, pch=16, cex=.8, ylim=c(-85, 20), ylab=expression(beta ~ 1), 
       xlab="Bandwidth Size", main="Democratic Primaries") #Local Linear RDDs
abline(h=0, lty=2)
dev.off()

rm(cf_data, Coef, Upper, Lower, SE, VarCovar, i, Local, num, Reg, Robust, X)



#######################################################################################################################
#######################################################################################################################
#######################################################################################################################
#######################################################################################################################
############################### Supplementary Materials
#######################################################################################################################
#######################################################################################################################

###################################
#### Balance Plots - Code adapted from Caughey and Sekhon (2011) Replication File
####################################
## Pre-treatment covariates
covars <- matrix(c('DWinNxt', 'Dem Win t + 1',
                   'DPctNxt', 'Dem % t + 1',
                   'DifDPNxt', 'Dem % Margin t + 1', 
                   'CQRating3', 'CQ Rating {-1, 0, 1}',
                   'DWinPrv', 'Dem Win t - 1',
                   'DPctPrv', 'Dem % t - 1',
                   'DifDPPrv', 'Dem % Margin t - 1',
                   'IncDWNOM1', 'Inc\'s D1 NOMINATE',
                   'DemInc', 'Dem Inc in Race',
                   'NonDInc', 'Rep Inc in Race',
                   'RExpAdv', 'Rep Experience Adv',
                   'DExpAdv', 'Dem Experience Adv',
                   'GovDem', 'Dem Governor',
                   'OpenSeat', 'Open Seat',
                   'VtTotPct', 'Voter Turnout %',
                   'GovWkPct', 'Pct Gov\'t Worker',
                   'UrbanPct', 'Pct Urban',
                   'BlackPct', 'Pct Black',
                   'ForgnPct', 'Pct Foreign Born',
                   'ElcSwing', 'Partisan Swing', 
                   'DSpndPct', 'Dem Spending %',
                   'DDonaPct', 'Dem Donation %',
                   'SoSDem', 'Dem Sec of State',
                   'DifPVDec', 'Dem Pres % Margin',
                   'DemOpen', 'Dem-held Open Seat',
                   'NonDOpen', 'Rep-held Open Seat', 
                   'PrvTrmsD', 'Dem\'s # Prev Terms',
                   'PrvTrmsO', 'Rep\'s # Prev Terms',
                   'IncumbentOppo', "Vs. Incumbent (Gen)",
                   'Primary_Oppo_Inc', "Vs. Incumbent (P)",
                   'redist', "Redistricting",
                   'year', "Year", 
                   'candnumber', "# of Primary Candidates",
                   'prez', "Pres. % T-1",
                   'PrezYear', "Presidential Elec Year",
                   "votep", "Party Vote T-1") , ncol = 2, byrow = TRUE)

########## Balance plot function
###### Function sends pdf to working directory --- 5% BW = -2.5% --> 2.5%
bp_func <- function(bandwidth, file_name, input_data, covs){
  ## Parameters for plotting
  r <- 10000 ; varline <- 6 ; nline <- 4 ; tline <- 2 ; cline <- 0
  data <- input_data
  bw <- abs(data$PrimMargin) <= bandwidth/2
  pdf(file_name , width = 7, height = 8)
  par(mar = c(2, 12, 1, 0))
  plot(x = NULL, y = NULL, xlim = c(0, 1), ylim = c(1, nrow(covs)), ylab = '', xlab = '', xaxt = "n", yaxt = "n", bty = 'n')
  mtext(text = c('Variable\nName', 'Valid\nCases', 'Treated\nMean', 'Control\nMean'), side = 2, font = 2, line = c(varline + 3, nline + .7, tline + .7, cline + 0.3), adj = .5, las = 2, at = nrow(covs) +1, cex = .7)
  ## For each covariate...
  for(i in 1:nrow(covs)) {
    print(covs[i, 2])
    print(sum(is.na(data[bw, covs[i, 1]])))
    aty <- nrow(covs) - i + 1
    print(aty)
    mtext(text = covs[i, 2], side = 2, line = varline, adj = 1, las = 2, at = aty, cex = .7)
    ## Number of valid cases
    nonna <- sum(!is.na(data[bw, covs[i, 1]]))
    ## Mean of treated
    meanT <- signif(mean(data[(bw & data$primary_winner == 1), covs[i, 1]], na.rm = TRUE), digits = 2)
    ## Adding/subtracting digits for presentation purposes
    if(abs(meanT) < 0.1) { meanT <- signif(meanT, digits = 1) }
    if(meanT %% 1 == 0 & abs(meanT) < 10) { meanT <- paste(meanT, '.0', sep = '') }
    if(abs(as.numeric(meanT)) >= .1 & abs(as.numeric(meanT)) < 1 & nchar(meanT) == 3) {meanT <- paste(meanT, '0', sep = '')}
    ## Mean of control
    meanC <- signif(mean(data[(bw & data$primary_winner == 0), covs[i, 1]],na.rm = TRUE), digits = 2)
    ## Presentation adjustments
    if(abs(meanC) < 0.1) {meanC <- signif(meanC, digits = 1)}
    if(meanC %% 1 == 0 & abs(meanC) < 10) { meanC <- paste(meanC, '.0', sep = '')}
    if(as.numeric(meanC) %% .1 == 0 & abs(as.numeric(meanC)) < 1) {meanC <- paste(meanC, '0', sep = '')}
    mtext(text = c(nonna, meanT, meanC), side = 2,line = c(nline, tline, cline),adj = 1, las = 2, at = aty, cex = .7) 
    ## Gray bands
    if(aty %% 2 == 1) { polygon(x = c(0, 0, 1, 1), y = c(aty - .5, aty + .5, aty + .5, aty - .5), border = FALSE,col = 'lightgray')}
    ## If the variable has three or more levels
    if(length(levels(factor(data[, covs[i, 1]]))) >= 3) {
      p1 <- pvalue(wilcox_test(data[, covs[i, 1]][bw] ~ factor(data$primary_winner[bw]), distribution = 'exact'))
      print(p1)
      sym <- 18
    } else
      ## If the variable is dichotomous
    {
      p1 <- fisher.test(x = factor(data[, covs[i, 1]][bw]), y = factor(data$primary_winner[bw]))$p.value
      sym <- 20
    }
    ## Plot p-value
    points(pch = sym, x = p1, y = aty)
  }
  segments(x0 = 0, x1 = 0, y0 = .5, y1 = nrow(covs) + .4)
  segments(x0 = 0, x1 = 1, y0 = .49, y1 = .49)
  segments(x0 = c(.05, .1), x1 = c(.05, .1), y0 = .5, y1 = nrow(covs) + .4, lty = 'dotted')
  mtext(side = 1, at = c(0, .05, .1, 1), text = c('0', '.05', '.1', '1'),cex = .7, line = -.75)
  mtext(side = 1, at = .5, text = 'p-value')
  #  mtext(side=3, at=.25, text=paste(pretitle, "Covariate Balance for ", bandwidth*100, "% Bandwidth", sep=""))
  dev.off()
}

#################################
###### Supplementary Figures 1 - 4 - Balance Plots
####################################

## Democrats
bp_func(.02, paste0(figure_dir,"Dem_Bal2.pdf"), dat[dat$party == "D",], covars[c(1:4, 6:8, 10:24, 26, 28:36),])
# bp_func(.03, paste0(figure_dir,"Dem_Bal3.pdf"), dat[dat$party == "D",], covars[c(1:8, 10:26, 28:36),])
# bp_func(.10, paste0(figure_dir,"Dem_Bal10.pdf"), dat[dat$party == "D",], covars[c(1:8, 10:26, 28:36),])

## Republicans
bp_func(.02, paste0(figure_dir,"Rep_Bal2.pdf"), dat[dat$party == "R",], covars[c(1:9, 11:24, 26:27, 29, 31:36),])
bp_func(.03, paste0(figure_dir,"Rep_Bal3.pdf"), dat[dat$party == "R",], covars[c(1:9, 11:24, 26:27, 29, 31:36),])
bp_func(.10, paste0(figure_dir,"Rep_Bal10.pdf"), dat[dat$party == "R",], covars[c(1:9, 11:24, 26:27, 29, 31:36),])
# *** Note: Dropping some vars as a result of insufficient variation within bandwidth (which means no observed difference)
dev.off()

############################################################################
############# Local Linear models of covariates
#####################################################################

#################################
###### Supplementary Figure 5 - Democrats Only
####################################
pdf(paste0(figure_dir, "Dem_Cov_LL.pdf"), height=6, width=16)
Sequence <- seq(.10, .25, .01)
par(mfrow=c(2, 4))

Form <- formula(inc ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + Y82to90 + Y92to00 + Y02to10)
Inc <- LocalLin(Seq=Sequence, Sub=dat$party == "D", Formula=Form, PredProb=F)
plotCI(x=Inc$Seq, y=Inc$Est, ui= Inc$High, li= Inc$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Incumbency", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

Form <- formula(prez ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + PrezYear + Y82to90 + Y92to00 + Y02to10)
Prez <- LocalLin(Seq=Sequence, Sub=dat$party == "D", Formula=Form, PredProb=F)
plotCI(x=Prez$Seq, y=Prez$Est, ui= Prez$High, li= Prez$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Presidential Vote Share", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

Form <- formula(candnumber ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + PrezYear  + Y82to90 + Y92to00 + Y02to10)
NumC <- LocalLin(Seq=Sequence, Sub=dat$party == "D", Formula=Form, PredProb=F)
plotCI(x=NumC$Seq, y=NumC$Est, ui= NumC$High, li= NumC$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Number of Candidates", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

#### Two-Party Vote Share --- Lots of Missing Data
Form <- formula(votep ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + PrezYear + Y82to90 + Y92to00 + Y02to10)
VoteP <- LocalLin(Seq=Sequence, Sub=dat$party == "D", Formula=Form, PredProb=F)
plotCI(x=VoteP$Seq, y=VoteP$Est, ui= VoteP$High, li= VoteP$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="District Vote Share (T-1)", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

###### Year
Form <- formula(year ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + PrezYear + Y82to90 + Y92to00 + Y02to10)
Year <- LocalLin(Seq=Sequence, Sub=dat$party == "D", Formula=Form, PredProb=F)
plotCI(x=Year$Seq, y=Year$Est, ui= Year$High, li= Year$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Year", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

#Redist
Form <- formula(redist ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + PrezYear + Y82to90 + Y92to00 + Y02to10)
redist <- LocalLin(Seq=Sequence, Sub=dat$party == "D", Formula=Form, PredProb=F)
plotCI(x=redist$Seq, y=redist$Est, ui= redist$High, li= redist$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Redistricting Year", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

#PrezYear
Form <- formula(PrezYear ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + Y82to90 + Y92to00 + Y02to10)
PYear <- LocalLin(Seq=Sequence, Sub=dat$party == "D", Formula=Form, PredProb=F)
plotCI(x=PYear$Seq, y=PYear$Est, ui= PYear$High, li= PYear$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Presidential Election Year", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

#IncumbentOppo
Form <- formula(IncumbentOppo ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + Y82to90 + Y92to00 + Y02to10)
IncOpp <- LocalLin(Seq=Sequence, Sub=dat$party == "D", Formula=Form, PredProb=F)
plotCI(x=IncOpp$Seq, y=IncOpp$Est, ui= IncOpp$High, li= IncOpp$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Incumbent Opponent in General", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)
dev.off()

#################################
###### Supplementary Figure 6 - Republicans Only
####################################

pdf(paste0(figure_dir, "Rep_Cov_LL.pdf"), height=6, width=16)
Sequence <- seq(.10, .25, .01)
par(mfrow=c(2, 4))

Form <- formula(inc ~ primary_winner + PrimMargin + primary_winner*PrimMargin + prez + PrezYear + Y82to90 + Y92to00 + Y02to10)
Inc <- LocalLin(Seq=Sequence, Sub=dat$party == "R", Formula=Form, PredProb=F)
plotCI(x=Inc$Seq, y=Inc$Est, ui= Inc$High, li= Inc$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Incumbency", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

Form <- formula(prez ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + PrezYear + Y82to90 + Y92to00 + Y02to10)
Prez <- LocalLin(Seq=Sequence, Sub=dat$party == "R", Formula=Form, PredProb=F)
plotCI(x=Prez$Seq, y=Prez$Est, ui= Prez$High, li= Prez$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Presidential Vote Share", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

Form <- formula(candnumber ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + PrezYear + Y82to90 + Y92to00 + Y02to10)
NumC <- LocalLin(Seq=Sequence, Sub=dat$party == "R", Formula=Form, PredProb=F)
plotCI(x=NumC$Seq, y=NumC$Est, ui= NumC$High, li= NumC$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Number of Candidates", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

#### Two Party Vote Share (NAs)
Form <- formula(votep ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + PrezYear + Y82to90 + Y92to00 + Y02to10)
VoteP <- LocalLin(Seq=Sequence, Sub=dat$party == "R", Formula=Form, PredProb=F)
plotCI(x=VoteP$Seq, y=VoteP$Est, ui= VoteP$High, li= VoteP$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="District Vote Share (T-1)", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

### Year
Form <- formula(year ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + PrezYear + Y82to90 + Y92to00 + Y02to10)
Year <- LocalLin(Seq=Sequence, Sub=dat$party == "R", Formula=Form, PredProb=F)
plotCI(x=Year$Seq, y=Year$Est, ui= Year$High, li= Year$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Year", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

#Redist
Form <- formula(redist ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + PrezYear + Y82to90 + Y92to00 + Y02to10)
redist <- LocalLin(Seq=Sequence, Sub=dat$party == "R", Formula=Form, PredProb=F)
plotCI(x=redist$Seq, y=redist$Est, ui= redist$High, li= redist$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Redistricting Year", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

#PrezYear
Form <- formula(PrezYear ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + Y82to90 + Y92to00 + Y02to10)
PYear <- LocalLin(Seq=Sequence, Sub=dat$party == "R", Formula=Form, PredProb=F)
plotCI(x=PYear$Seq, y=PYear$Est, ui= PYear$High, li= PYear$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Presidential Election Year", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)

#IncumbentOppo
Form <- formula(IncumbentOppo ~ primary_winner + PrimMargin + primary_winner*PrimMargin + inc + prez + PrezYear + Y82to90 + Y92to00 + Y02to10)
IncOpp <- LocalLin(Seq=Sequence, Sub=dat$party == "R", Formula=Form, PredProb=F)
plotCI(x=IncOpp$Seq, y=IncOpp$Est, ui= IncOpp$High, li= IncOpp$Low, pch=16, cex=.8, col=alpha("black", .75),
       main="Incumbent Opponent in General", xlab="Bandwidth", ylab=expression(beta ~ 1)); abline(h=0, lty=2)
dev.off()

#################################
###### Supplementary Figure 7 - Histograms by Party and Bandwidth
####################################

pdf(paste0(figure_dir, "FemaleMargins_byParty.pdf"), height=8.4, width=14)
par(mfrow=c(2, 3), mar=c(5.1, 4.1, 3.4, 2.1))
hist(x = dat[abs(dat$PrimMargin) <= .1 & dat$party == "R",]$PrimMargin, breaks = seq(-.1, .1, 0.025),
     freq = TRUE, labels = TRUE, right = FALSE, xlim = c(-.1, .1), col = alpha("grey", 1), main = "", ylim=c(0, 25),
     xlab = "Female Candidate Primary Margin (%)", ylab = 'Frequency Count in 2.5% Bins', las = 1)
abline(v = 0, lwd = 2, lty=2) 
hist(x = dat[abs(dat$PrimMargin) <= .1  & dat$party == "R",]$PrimMargin, breaks = seq(-.1, .1, 0.01),
     freq = TRUE, labels = TRUE, right = FALSE, xlim = c(-.1, .1), col = alpha("grey", 1), main = "Republican Candidates", 
     ylim=c(0, 12),xlab = "Female Candidate Primary Margin (%)", ylab = 'Frequency Count in 1.0% Bins', las = 1)
abline(v = 0, lwd = 2, lty=2) 
hist(x =dat[abs(dat$PrimMargin) <= .1  & dat$party == "R",]$PrimMargin, breaks = seq(-.1, .1, 0.005),
     freq = TRUE, labels = TRUE, right = FALSE, xlim = c(-.1, .1), col = alpha("grey", 1), main = "", ylim=c(0, 8),
     xlab = "Female Candidate Primary Margin (%)", ylab = 'Frequency Count in 0.5% Bins', las = 1)
abline(v = 0, lwd = 2, lty=2)
hist(x = dat[abs(dat$PrimMargin) <= .1 & dat$party == "D",]$PrimMargin, breaks = seq(-.1, .1, 0.025),
     freq = TRUE, labels = TRUE, right = FALSE, xlim = c(-.1, .1), col = alpha("grey", .4), main = "", ylim=c(0, 25),
     xlab = "Female Candidate Primary Margin (%)", ylab = 'Frequency Count in 2.5% Bins', las = 1)
abline(v = 0, lwd = 2, lty=2) 
hist(x = dat[abs(dat$PrimMargin) <= .1  & dat$party == "D",]$PrimMargin, breaks = seq(-.1, .1, 0.01),
     freq = TRUE, labels = TRUE, right = FALSE, xlim = c(-.1, .1), col = alpha("grey", .4), main = "Democratic Candidates",
     ylim=c(0, 12), xlab = "Female Candidate Primary Margin (%)", ylab = 'Frequency Count in 1.0% Bins', las = 1)
abline(v = 0, lwd = 2, lty=2) 
hist(x =dat[abs(dat$PrimMargin) <= .1  & dat$party == "D",]$PrimMargin, breaks = seq(-.1, .1, 0.005),
     freq = TRUE, labels = TRUE, right = FALSE, xlim = c(-.1, .1), col = alpha("grey", .4), main = "", ylim=c(0, 8),
     xlab = "Female Candidate Primary Margin (%)", ylab = 'Frequency Count in 0.5% Bins', las = 1)
abline(v = 0, lwd = 2, lty=2) 
dev.off()

#################################
###### Supplementary Figure 8 - McCrary Sorting Tests
####################################

pdf(paste0(figure_dir,"MC_Sorting.pdf"), height=4, width=8)
par(mfrow=c(1,2), cex.axis = 1.5, cex.lab = 1.5)
DCdensity(dat[dat$party == "D",]$PrimMargin, 0, verbose=T); title(main = 'Democratic Primaries')
DCdensity(dat[dat$party == "R",]$PrimMargin, 0, verbose=T); title(main = 'Republican Primaries')
dev.off()


################################################################################
#### Supplementary Figure 9: Effect Heterogeneity by District Partisanship (Republicans)
################################################################################
# ** Referencened in Footnote 17

rep_interact_win <- lm(GenWin ~ primary_winner*prez + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner, data=dat[dat$party == "R",])
rep_interact_win2 <- lm(GenWin ~ primary_winner*prez + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner, data=dat[dat$party == "R" & abs(dat$PrimMargin) < .2,])
rep_interact_win3 <- lm(GenWin ~ primary_winner*prez + PrimMargin + primary_winner*PrimMargin + I(PrimMargin^2)*primary_winner + PrezYear + inc + Y82to90 + Y92to00 + Y02to10, data=dat[dat$party == "R",])

base <- interplot::interplot(rep_interact_win, "primary_winner", "prez", hist = TRUE) + geom_hline(aes(yintercept = 0), lty = 2, lwd = .5)
sub <- interplot::interplot(rep_interact_win2, "primary_winner", "prez", hist = TRUE) + geom_hline(aes(yintercept = 0), lty = 2, lwd = .5)
controls <- interplot::interplot(rep_interact_win3, "primary_winner", "prez", hist = TRUE) + geom_hline(aes(yintercept = 0), lty = 2, lwd = .5)

base <- base + xlab("Presidential Vote Share, T-1") + ylab("General Election Win Probability") + ggtitle("No Controls")
sub <- sub + xlab("Presidential Vote Share, T-1") + ylab("General Election Win Probability") + ggtitle("No Controls - Close races")
controls <- controls + xlab("Presidential Vote Share, T-1") + ylab("General Election Win Probability") + ggtitle("Full Controls")

pdf(paste0(figure_dir, "Prez_Interact.pdf"), height = 4, width = 12)
grid.arrange(base, sub, controls, 
             nrow = 1, 
             top = "Effect of Nominating a Female Candidate on General Election Win Probability Across Levels of District Partisanship")
dev.off()

################################################################################
#### Supplementary Table 1: Effect Heterogeneity by District Partisanship
################################################################################

stargazer(rep_interact_win, rep_interact_win2, rep_interact_win3,
          omit.stat=c("adj.rsq", "f", "ser"), style="apsr",  #type = "text",
          title="Effect of Nominating a Female Candidate on General Election Win Probability Across Levels of District Partisanship",
          covariate.labels=c("Female Win", "Pres. Vote T-1",  "Margin", "Margin^2","Pres. Election Year",
                             "Incumbent","1982 - 1990", "1992 - 2000", "2000 - 2010","Female Win * Pres. Vote T-1",
                             "Female Win * Margin", "Female Win * Margin^2", "Constant"))


